#include "stegit_math.h"

#include <gsl/gsl_interp.h>

const unsigned int CStegitMath::sLogApproxRounds(100);

double CStegitMath::pow(double pBase, int pPower)
{
	double lResult(pBase);

	for (int i = 1; i < pPower; i++) {
		lResult *= pBase;
	}

	return lResult;
}

std::complex<double> CStegitMath::pow(double pBase, std::complex<double> pPower)
{
	// a^z = exp(z * log(a))
	//
	// exp(z * log(a)) = e^x * cos y + i * e^x * sin y
	// e^x = sinh(x) + cosh(x)

	double lLog(CStegitMath::log(pBase));
	std::complex<double> lExp(real(pPower) * lLog, imag(pPower) * lLog);

	double lEexp(cosh(real(lExp)) + sinh(real(lExp))),
		   lCos(cos(imag(lExp))),
		   lSin(sin(imag(lExp)));

	return std::complex<double>(lEexp * lCos, lEexp * lSin);
}

double CStegitMath::log(double pValue)
{
	double lResult(0.0), lDivisor, lNumerator, lPower;

	// approximation for logarithm function
	for (unsigned int i = 0; i < CStegitMath::sLogApproxRounds; i++) {
		lPower = (2*i) + 1;
		lDivisor = CStegitMath::pow((pValue-1), lPower);
		lNumerator = CStegitMath::pow((pValue+1), lPower);

		lResult += lDivisor / (lPower * lNumerator);
	}

	return (2*lResult);
}

// convert double to c++ complex
void CStegitMath::convert(double* pInput, std::complex<double>* pOutput,
		unsigned int pLength)
{
#if 0
	// this works for 8-bit pcm samples

	char lChar;
	double lDouble;

	for (unsigned int i = 0; i < pLength; i++) {
		lChar = pInput[i];
		lDouble = 0.0;

		if (lChar < 0) {
			lDouble = lChar + 128.0;
		} else if (lChar > 0) {
			lDouble = lChar - 128.0;
		}

		pOutput[i] = std::complex<double>(lDouble / 127.0, 0.0);
	}
#endif

	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i] = std::complex<double>(pInput[i], 0.0);
	}
}

// convert c++ complex to double
void CStegitMath::convert(std::complex<double>* pInput, double* pOutput,
		unsigned int pLength)
{
#if 0
	// this works for 8-bit pcm samples

	double lDouble;

	for (unsigned int i = 0; i < pLength; i++) {
		lDouble = real(pInput[i]) * 127.0;

		if (real(pInput[i]) < 0.0) {
			lDouble += 128.0;
		} else { // >= 0.0
			lDouble -= 128.0;
		}

		pOutput[i] = lDouble;
	}
#endif

	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i] = real(pInput[i]);
	}
}

// convert c++ complex to fftw complex
void CStegitMath::convert(std::complex<double>* pInput, fftw_complex* pOutput,
		unsigned int pLength)
{
	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i][0] = real(pInput[i]); // real part
		pOutput[i][1] = imag(pInput[i]); // imaginary part
	}
}

// convert double to fftw complex
void CStegitMath::convert(double* pInput, fftw_complex* pOutput,
		unsigned int pLength)
{
	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i][0] = pInput[i]; // real part
		pOutput[i][1] = 0.0; // imaginary part
	}
}

// convert fftw complex to c++ complex
void CStegitMath::convert(fftw_complex* pInput, std::complex<double>* pOutput,
		unsigned int pLength)
{
	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i] = std::complex<double>(pInput[i][0], pInput[i][1]);
	}
}

// convert fftw complex to double
void CStegitMath::convert(fftw_complex* pInput, double* pOutput,
		unsigned int pLength)
{
	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i] = pInput[i][0];
	}
}

// convert valarray<double> to fftw complex
void CStegitMath::convert(std::valarray<double>& pInput, fftw_complex* pOutput,
		unsigned int pLength)
{
	for (unsigned int i = 0; i < pLength; i++) {
		pOutput[i][0] = pInput[i]; // real part
		pOutput[i][1] = 0.0; // imaginary part
	}
}

bool CStegitMath::interpolate
(
		double* pOut,
		unsigned int pOutLen,
		double* pIn,
		unsigned int pInLen,
		double* pInTmp,
		EInterpolationType pType,
		bool pExtrapolation
)
{
	if(!pOut || !pIn || !pInTmp || pOutLen == 0 || pInLen == 0) {
		return false;
	}

	const gsl_interp_type* lInterpType;
	switch(pType) {
		case INTERP_TYPE_LINEAR:
			lInterpType = gsl_interp_linear;
			break;
		case INTERP_TYPE_POLYNOMIAL:
			lInterpType = gsl_interp_polynomial;
			break;
		case INTERP_TYPE_CSPLINE:
			lInterpType = gsl_interp_cspline;
			break;
		case INTERP_TYPE_CSPLINE_PERIODIC:
			lInterpType = gsl_interp_cspline_periodic;
			break;
		case INTERP_TYPE_AKIMA:
			lInterpType = gsl_interp_akima;
			break;
		case INTERP_TYPE_AKIMA_PERIODIC:
			lInterpType = gsl_interp_akima_periodic;
			break;
		default:
			return false;
	}

	double lXSpace;
	if(pExtrapolation) {
		lXSpace = 1.0 / static_cast<double>(pInLen);
	}
	else {
		lXSpace = 1.0 / static_cast<double>(pInLen - 1);
	}
	double lXOffset = 0.0;
	for(unsigned int i = 0; i < pInLen; i++) {
		pInTmp[i] = lXOffset;
		lXOffset += lXSpace;
	}

	// allocate the gsl objects
	gsl_interp_accel* lAccel = gsl_interp_accel_alloc();
	gsl_interp* lInterp = gsl_interp_alloc(lInterpType, pInLen);

	// inialise the gsl objects
	gsl_interp_init(lInterp, pInTmp, pIn, pInLen);

	if(pExtrapolation) {
		// the value of lXOffset is the x-length of input samples ( < 1.0 )
		lXSpace = lXOffset / static_cast<double>(pOutLen);
	}
	else {
		lXSpace = 1.0 / static_cast<double>(pOutLen - 1);
	}
	lXOffset = 0.0;
	for(unsigned int i = 0; i < pOutLen; i++) {
		pOut[i] = gsl_interp_eval(lInterp, pInTmp, pIn, lXOffset, lAccel);
		lXOffset += lXSpace;
	}
	gsl_interp_free(lInterp);
	gsl_interp_accel_free(lAccel);

	return true;
}

bool CStegitMath::interpolate
(
		double* pOut,
		unsigned int pOutLen,
		double* pIn,
		unsigned int pInLen,
		EInterpolationType pType,
		bool pExtrapolation
)
{
	if(pInLen == 0) {
		return false;
	}
	double* lInTmp = new double [pInLen];
	bool rv(interpolate(pOut, pOutLen, pIn, pInLen, lInTmp, pType, pExtrapolation));
	delete [] lInTmp;
	return rv;
}

bool CStegitMath::interpolate
(
		double* pOut,
		unsigned int pOutLen,
		double pOutLenTime,
		double pOutOffsetTime,
		double* pIn,
		unsigned int pInLen,
		double pInLenTime,
		EInterpolationType pType
)
{
	if(!pOut || !pIn || pOutLen == 0 || pInLen == 0 || pOutLenTime <= 0.0 ||
			pInLenTime <= 0.0) {
		return false;
	}

	const gsl_interp_type* lInterpType;
	switch(pType) {
		case INTERP_TYPE_LINEAR:
			lInterpType = gsl_interp_linear;
			break;
		case INTERP_TYPE_POLYNOMIAL:
			lInterpType = gsl_interp_polynomial;
			break;
		case INTERP_TYPE_CSPLINE:
			lInterpType = gsl_interp_cspline;
			break;
		case INTERP_TYPE_CSPLINE_PERIODIC:
			lInterpType = gsl_interp_cspline_periodic;
			break;
		case INTERP_TYPE_AKIMA:
			lInterpType = gsl_interp_akima;
			break;
		case INTERP_TYPE_AKIMA_PERIODIC:
			lInterpType = gsl_interp_akima_periodic;
			break;
		default:
			return false;
	}

	double lXSpace(pInLenTime / static_cast<double>(pInLen - 1));

	double lXOffset = 0.0;
	double* lXIn = new double [pInLen];
	for(unsigned int i = 0; i < pInLen; i++) {
		lXIn[i] = lXOffset;
		lXOffset += lXSpace;
	}

	// allocate the gsl objects
	gsl_interp_accel* lAccel = gsl_interp_accel_alloc();
	gsl_interp* lInterp = gsl_interp_alloc(lInterpType, pInLen);

	// inialise the gsl objects
	gsl_interp_init(lInterp, lXIn, pIn, pInLen);

	lXSpace = pOutLenTime / static_cast<double>(pOutLen - 1);
	lXOffset = pOutOffsetTime;
	for(unsigned int i = 0; i < pOutLen; i++) {
		pOut[i] = gsl_interp_eval(lInterp, lXIn, pIn, lXOffset, lAccel);
		lXOffset += lXSpace;
	}
	gsl_interp_free(lInterp);
	gsl_interp_accel_free(lAccel);

	delete [] lXIn;
	return true;
}

/**
 * Create a 16 bit CRC
 * in_cksum --
 *	Checksum routine for Internet Protocol family headers (C Version)
 * http://www.cs.utah.edu/~swalton/listings/sockets/programs/part4/chap18/ping.c
 */
uint16_t
CStegitMath::inCksum(void *pAddr, int pLen)
{
	register int nleft = pLen;
	register uint16_t *w = static_cast<uint16_t*>(pAddr);
	register int sum = 0;
	uint16_t answer = 0;

	/*
	 * Our algorithm is simple, using a 32 bit accumulator (sum), we add
	 * sequential 16 bit words to it, and at the end, fold back all the
	 * carry bits from the top 16 bits into the lower 16 bits.
	 */
	while (nleft > 1)  {
		sum += *w++;
		nleft -= 2;
	}

	/* mop up an odd byte, if necessary */
	if (nleft == 1) {
		*(unsigned char *)(&answer) = *(unsigned char *)w ;
		sum += answer;
	}

	/* add back carry outs from top 16 bits to low 16 bits */
	sum = (sum >> 16) + (sum & 0xffff);	/* add hi 16 to low 16 */
	sum += (sum >> 16);			/* add carry */
	answer = ~sum;				/* truncate to 16 bits */
	return answer;
}

// ================================ CRC-8 ==========================
/*
see:
http://www.koders.com/cpp/fid9C544B36B8C41721691790197D38DAC91D2C29EF.aspx?s=md5
*/

#define CRC8_INIT_VALUE		0x0000
#define CRC8_XOR_VALUE		0x0000

//#ifdef CREATE_CRC_TABLE
#if 0

static byte sCRCTable[256];

/*
   Generate a table for a byte-wise 8-bit CRC calculation on the polynomial:
   x^8 + x^2 + x^1 + x^0
*/

void CStegitMath::make_crc_table( void ) {
	int i, j;
	unsigned long poly, c;
	/* terms of polynomial defining this crc (except x^8): */
	static const byte p[] = {0,1,2};

	/* make exclusive-or pattern from polynomial (0x07) */
	poly = 0L;
	for ( i = 0; i < sizeof( p ) / sizeof( byte ); i++ ) {
		poly |= 1L << p[i];
	}

	for ( i = 0; i < 256; i++ ) {
		c = i;
		for ( j = 0; j < 8; j++ ) {
			c = ( c & 0x80 ) ? poly ^ ( c << 1 ) : ( c << 1 );
		}
		sCRCTable[i] = (byte) c;
	}
}

#else

/*
  Table of CRC-8's of all single-byte values (made by make_crc_table)
*/
unsigned char CStegitMath::sCRCTable[256] = {
	0x00, 0x07, 0x0E, 0x09, 0x1C, 0x1B, 0x12, 0x15,
	0x38, 0x3F, 0x36, 0x31, 0x24, 0x23, 0x2A, 0x2D,
	0x70, 0x77, 0x7E, 0x79, 0x6C, 0x6B, 0x62, 0x65,
	0x48, 0x4F, 0x46, 0x41, 0x54, 0x53, 0x5A, 0x5D,
	0xE0, 0xE7, 0xEE, 0xE9, 0xFC, 0xFB, 0xF2, 0xF5,
	0xD8, 0xDF, 0xD6, 0xD1, 0xC4, 0xC3, 0xCA, 0xCD,
	0x90, 0x97, 0x9E, 0x99, 0x8C, 0x8B, 0x82, 0x85,
	0xA8, 0xAF, 0xA6, 0xA1, 0xB4, 0xB3, 0xBA, 0xBD,
	0xC7, 0xC0, 0xC9, 0xCE, 0xDB, 0xDC, 0xD5, 0xD2,
	0xFF, 0xF8, 0xF1, 0xF6, 0xE3, 0xE4, 0xED, 0xEA,
	0xB7, 0xB0, 0xB9, 0xBE, 0xAB, 0xAC, 0xA5, 0xA2,
	0x8F, 0x88, 0x81, 0x86, 0x93, 0x94, 0x9D, 0x9A,
	0x27, 0x20, 0x29, 0x2E, 0x3B, 0x3C, 0x35, 0x32,
	0x1F, 0x18, 0x11, 0x16, 0x03, 0x04, 0x0D, 0x0A,
	0x57, 0x50, 0x59, 0x5E, 0x4B, 0x4C, 0x45, 0x42,
	0x6F, 0x68, 0x61, 0x66, 0x73, 0x74, 0x7D, 0x7A,
	0x89, 0x8E, 0x87, 0x80, 0x95, 0x92, 0x9B, 0x9C,
	0xB1, 0xB6, 0xBF, 0xB8, 0xAD, 0xAA, 0xA3, 0xA4,
	0xF9, 0xFE, 0xF7, 0xF0, 0xE5, 0xE2, 0xEB, 0xEC,
	0xC1, 0xC6, 0xCF, 0xC8, 0xDD, 0xDA, 0xD3, 0xD4,
	0x69, 0x6E, 0x67, 0x60, 0x75, 0x72, 0x7B, 0x7C,
	0x51, 0x56, 0x5F, 0x58, 0x4D, 0x4A, 0x43, 0x44,
	0x19, 0x1E, 0x17, 0x10, 0x05, 0x02, 0x0B, 0x0C,
	0x21, 0x26, 0x2F, 0x28, 0x3D, 0x3A, 0x33, 0x34,
	0x4E, 0x49, 0x40, 0x47, 0x52, 0x55, 0x5C, 0x5B,
	0x76, 0x71, 0x78, 0x7F, 0x6A, 0x6D, 0x64, 0x63,
	0x3E, 0x39, 0x30, 0x37, 0x22, 0x25, 0x2C, 0x2B,
	0x06, 0x01, 0x08, 0x0F, 0x1A, 0x1D, 0x14, 0x13,
	0xAE, 0xA9, 0xA0, 0xA7, 0xB2, 0xB5, 0xBC, 0xBB,
	0x96, 0x91, 0x98, 0x9F, 0x8A, 0x8D, 0x84, 0x83,
	0xDE, 0xD9, 0xD0, 0xD7, 0xC2, 0xC5, 0xCC, 0xCB,
	0xE6, 0xE1, 0xE8, 0xEF, 0xFA, 0xFD, 0xF4, 0xF3
};

#endif

void CStegitMath::CRC8_InitChecksum(unsigned char &crcvalue)
{
	crcvalue = CRC8_INIT_VALUE;
}

void CStegitMath::CRC8_Update(unsigned char &crcvalue, const unsigned char data)
{
	crcvalue = sCRCTable[crcvalue ^ data];
}

void CStegitMath::CRC8_UpdateChecksum(unsigned char &crcvalue, const void *data, int length)
{
	unsigned char crc;
	const unsigned char *buf = (const unsigned char *) data;

	crc = crcvalue;
	while( length-- ) {
		crc = sCRCTable[crc ^ *buf++];
	}
	crcvalue = crc;
}

void CStegitMath::CRC8_FinishChecksum(unsigned char &crcvalue)
{
	crcvalue ^= CRC8_XOR_VALUE;
}

unsigned char CStegitMath::CRC8_BlockChecksum(const void *data, int length)
{
	unsigned char crc;

	CRC8_InitChecksum(crc);
	CRC8_UpdateChecksum(crc, data, length);
	CRC8_FinishChecksum(crc);
	return crc;
}
